{
"cells": [
{
"cell_type": "markdown",
"id": "8f039ec7",
"metadata": {},
"source": [
"# Differentiation\n",
"**강좌**: *수치해석*"
]
},
{
"cell_type": "markdown",
"id": "13f89f04",
"metadata": {},
"source": [
"## Finite Difference\n",
"Taylor expansion 을 이용하면 수치 미분을 구할 수 있다.\n",
"\n",
"### First-order derivative\n",
"\n",
":::{figure-md} fd\n",
"
\n",
"\n",
"Finite Difference (from Wikipedia)\n",
"::: \n",
"\n",
"- Forward Difference\n",
"\n",
"$$\n",
"f'(x_j) = \\frac{f(x_{j+1}) - f(x_j)}{\\Delta x} + O(\\Delta x) \\approx \\frac{\\Delta f(x)}{\\Delta x}\n",
"$$\n",
"\n",
"- Backward Differnce\n",
"\n",
"$$\n",
"f'(x_j) = \\frac{f(x_{j}) - f(x_{j-1})}{\\Delta x} + O(\\Delta x) \\approx \\frac{\\nabla f(x)}{\\Delta x}\n",
"$$\n",
"\n",
"- Central Difference\n",
" - 다음 두 식을 빼보자.\n",
" \n",
"$$\n",
"f(x_{j+1}) = f(x_j) + \\Delta x f'(x_j) + \\frac{(\\Delta x)^2}{2} f''(x_j) + \\frac{(\\Delta x)^3}{3!} f'''(x_j) + O((\\Delta x)^2)\n",
"$$\n",
"$$\n",
"f(x_{j-1}) = f(x_j) - \\Delta x f'(x_j) + \\frac{(\\Delta x)^2}{2} f''(x_j) - \\frac{(\\Delta x)^3}{3!} f'''(x_j) + O((\\Delta x)^3)\n",
"$$\n",
"\n",
" - 그 결과는 다음과 같다.\n",
" \n",
"$$\n",
"f'(x_j) = \\frac{f(x_{j+1}) - f(x_{j-1})}{ 2\\Delta x} + O((\\Delta x)^2) \\approx \\frac{\\delta f(x)}{2 \\Delta x}\n",
"$$\n",
"\n",
"### Second-order derivative\n",
"- Central Difference \n",
" - 위의 두 식을 더한 후 $2 f(x_j)$ 를 빼보자.\n",
" \n",
"$$\n",
"f''(x_j) = \\frac{f(x_{j+1}) + 2 f(x_j) - f(x_{j-1})}{ (\\Delta x)^2} + O((\\Delta x)^2)\n",
"$$ \n",
"\n",
"#### 제차분식\n",
"- Second-order Central\n",
"\n",
" $$\n",
" f''(x) \\approx \\frac{\\delta^2 f(x) }{ (\\Delta x)^2} = \\frac{\\Delta f(x) - \\nabla f(x) }{\\Delta x} = \\frac{f(x_{j+1}) + 2 f(x_j) - f(x_{j-1})}{ (\\Delta x)^2} \n",
" $$\n",
" \n",
" * $E=O((\\Delta x)^2)$\n",
"\n",
"- Second-order Forward\n",
"\n",
" $$\n",
" f''(x) \\approx \\frac{\\Delta^2 f(x) }{ (\\Delta x)^2} = \\frac{\\Delta f(x+\\Delta x) - \\Delta f(x) }{\\Delta x} = \\frac{f(x_{j+2}) - 2 f(x_{j+1}) + f(x_{j})}{ (\\Delta x)^2}\n",
" $$\n",
" \n",
" * $E=O((\\Delta x))$\n",
"\n",
"- Second-order Backward\n",
"\n",
" $$\n",
" f''(x) \\approx \\frac{\\nabla^2 f(x) }{ (\\Delta x)^2} = \\frac{\\nabla f(x) - \\nabla f(x - \\Delta x) }{\\Delta x} = \\frac{f(x_{j}) - 2 f(x_{j-1}) + f(x_{j-2})}{ (\\Delta x)^2}\n",
" $$\n",
" \n",
" * $E=O((\\Delta x))$\n",
"\n",
"### General approaches\n",
"- 일반적인 방법\n",
" - 여러 Taylor expansion의 합차를 이용해서 원하는 미분항의 근사식을 구한다.\n",
" \n",
"\n",
"### 예제\n",
"$f(x)=\\sin(x)$ 에 대해서 수치 미분 결과를 비교하자."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "8929c44f",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"from matplotlib import pyplot as plt\n",
"\n",
"import numpy as np\n",
"\n",
"plt.style.use('ggplot')\n",
"plt.rcParams['figure.dpi'] = 150"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "2f8012f3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact: 0.8090169943749475\n",
"0.7996517731793618 0.8181160727815516 0.8088839229804567\n"
]
}
],
"source": [
"x = 0.2*np.pi\n",
"dx = 0.01*np.pi\n",
"f = np.sin\n",
"\n",
"# First-order forward \n",
"fdf = (f(x+dx) - f(x)) / dx\n",
"\n",
"# First-order backward\n",
"bdf = (f(x) - f(x-dx)) / dx\n",
"\n",
"# First-order central\n",
"cdf = (f(x+dx) - f(x-dx)) / (2*dx)\n",
"\n",
"print(\"Exact: \", np.cos(x))\n",
"print(fdf, bdf, cdf)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "55849ca3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, '$\\\\epsilon$')"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAAKOCAYAAADu7ltwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAABcSAAAXEgFnn9JSAAB5HElEQVR4nO3deXxM9/7H8fcQSZgggqBpQxJLbHVr1wVXaRFUbWm1WqXVH1q0t3XvLapKl9tbbheqC9pLURWlJUppi9pprLU09qW2CMIQkWR+f+TONGOSyD5zZl7Px+M+bnO2+cx8ZSbv+ZzzPSar1WoVAAAAAMBwSri6AAAAAABA/hDoAAAAAMCgCHQAAAAAYFAEOgAAAAAwKAIdAAAAABgUgQ4AAAAADIpABwAAAAAGRaADAAAAAIMi0AEAAACAQRHoAAAAAMCgCHQAAAAAYFAEOgAAAAAwKAIdAAAAABgUgQ4AAAAADMrH1QWgYE6fPi2r1Vosj1WpUiVJUkJCQrE8HooPY+uZGFfPxdh6JsbVczG2ninzuJpMJlWtWtUldRDoDM5qtRZboMv8mPBMjK1nYlw9F2PrmRhXz8XYeiZXjyunXAIAAACAQRHoAAAAAMCgCHQAAAAAYFAEOgAAAAAwKAIdAAAAABgUgQ4AAAAADIrbFgAAAENy9VThhc32fDzteYGxNSqTyeTqEnKFQAcAAAwjPT1d165dU3Jyssf9cXzhwgVJGc8RnoWxNSaTySR/f3+VLl1aJUq474mNBDoAAGAIVqtVly5dko+Pj8qXL6+SJUu6uqRC5eOT8WdZamqqiytBYWNsjSktLU3Xrl3TpUuXFBgY6LYdOwIdAAAwhGvXrqlEiRIKCAhw2z+sCsL2nDzxuXk7xtaYfHx8FBAQoKSkJF27dk1lypRxdUlZct/eIQAAQCYpKSny9/fnj2IAxcZ22mVKSoqrS8kWgQ4AABhCamqqSpUq5eoyAHiZUqVKufXpsgQ6AADg9mwToNCdA1DcbO877joRE4EOAAAAAAyKQAcAAAAABkWgAwAAAACDItABAAAAgEER6AAAAAwqJCQkx//16tXL1SW6pZCQELVo0aJQjrV+/XqFhIRoxIgRDssnTpyokJAQzZs3z2mfnTt36tFHH1XdunXtY3X8+HFJ0tWrVzVmzBg1bdpUoaGhCgkJ0cSJEwulVngmbiwOAABgcL17985yec2aNYu5EtzKlStX9NRTT+nMmTNq1aqVbrvtNplMJpnNZknSW2+9pRkzZqhGjRrq2rWrSpUqpfr167u4argzAh0AAIDBvffee64uATd56qmn9NBDDyk4ONhh+fbt23X69Gn17NlTH3zwgdN+y5Ytk7+/v1asWKEyZcoUV7kwME65BAAAbsl6JUlpU9+SNfGcq0sB8iwoKEg1a9ZUuXLlHJafOnVKklS9evUs9zt16pQqVapEmEOuEegAAIDbsVouK33SGB38/ZjS3x0l63lCXWE4efKkRo4cqebNmyssLEx33nmnnn76aW3fvt1p2+PHj9uvw7t8+bLGjRunli1bqnr16nr11Vc1b948hYSEaNKkSQ77nT9/XrfffrtCQkL09ddfO6zbsWOHQkJCNGjQIPuyS5cuacaMGerbt6+9rvr16+uxxx7TmjVrsnwevXr1sl93tnDhQnXp0kW1a9dW3bp17dtcvXpVb7zxhpo1a6bw8HC1bt1an3zySb5vDn38+HENGTJE9evXV61atdStWzetXLky2+1vvobu2LFjDtfaTZo0yX793IgRI+zPyWq16sSJEw7XQmZ2/vx5vf7667rvvvsUHh6uevXq6fHHH9fGjRudash8fd/Zs2f10ksvqUmTJgoNDdVnn33m8NxGjhypFi1aKCwsTA0bNtQzzzyjPXv2OB3TNu4TJ07UyZMnNXToUDVs2FARERHq1KmTfvjhh2xfk99//10vvPCCfZwbNWqkHj16aNq0aU7bWiwW/ec//9H999+viIgI1alTRz179tSyZcuyPb634pTLIvbzzz9rzZo1OnbsmG7cuKFq1aqpS5cuuu+++1xdGgAAbslquaL0SWO0PDVYnzTtqejDPyh66y9S0zauLs3Q9u7dqz59+igxMVE1a9ZUp06ddPLkSX3//fdasWKFJk+erK5duzrtl5ycrJ49e+rkyZNq2bKlGjZsqMDAQLVq1UpSRmh48cUX7dtv2LDBHprWr1+vPn36OKyTZN9XkuLi4jRmzBjddtttCg8PV5MmTXTy5EmtXr1aq1ev1rvvvqtHHnkky+f04Ycfau7cuWrWrJnat2+vP/74Q5J0/fp1Pfroo9q6dauCgoLUvn17WSwWvfXWWzp69GieX7sjR47ooYceUkJCgsLDw9WwYUMdO3ZM/fv3V79+/XJ1DLPZrN69e+vIkSPasmWL6tWrZ782rnnz5qpZs6Zuv/12zZ8/X2XKlFFUVJTTMQ4cOKDo6GidPn1aNWrUULt27XThwgWtW7dOq1ev1gcffKCHH37Yab/z58+rc+fOSktLU7NmzXT9+nWVLl1akrR582Y98cQTunz5surUqaMOHTro9OnT+v777/XTTz9p5syZuueee5yOefz4cXXu3Fl+fn5q3ry5zp07p19//VUDBw7Ul19+qTZtHH9fFy9erOHDh+v69euqU6eOmjZtqosXL2r//v0aO3asnn76afu2586dU58+ffT777+ratWqat26ta5du2Y//j//+U8999xzuXrdvQGBrojt2rVLTZo00WOPPaaAgABt3rxZkydPVsmSJXX33Xe7ujwAANyK9eoVpf/nVX2fWkWf1cn4w3Re2AMyV68k5z8p/7eP1SpdsxRfkYWltFkmk6lYHspqter5559XYmKinnvuOf3jH/+wP/aSJUs0ePBgvfTSS2rZsqUqV67ssO+2bdvUpEkTrV+/XuXLl3dYFxISori4OCUnJ8vf319SRmgzmUyqVauWPcDZrF+/XpJjoIuIiNCiRYvUrFkzh213796tPn366LXXXlPXrl3tk4ZktmDBAn399dcOx5OkTz/9VFu3btVdd92lOXPm2E973LVrV7YTyOTklVdeUUJCgp588klNmDBBJUpknOQ2Z84cvfzyy7k6RsWKFfXee+9p3rx52rJlizp27Ki//e1vTtvNnz9fQUFBTtdFpqWl6dlnn9Xp06f1+uuva8CAAfYx3L17tx555BGNHDlS9913nypVquSw708//aROnTpp8uTJ9nGSpMuXL+vZZ59VcnKyPvnkE3Xp0sW+bs2aNXryySc1bNgwbdiwQb6+vk51DhgwQGPHjpWPT0akmDZtmsaOHav333/fIdAdOnRII0aMUHp6uqZOnapu3brZ16Wnp+vHH390OPaLL76o33//XUOGDNHIkSNVqlQpSdLRo0fVt29fvfPOO2rXrp3q1at3y9fdG3hUoDt06JB27typAwcOKD4+XhcuXFCpUqU0e/bsHPdLSUnRokWLtG7dOiUkJCggIECNGjVSdHS0KlasWKCahg0b5vBzt27d9Ntvv2n9+vUEOgAAMrFetSj9P2O1NDVY02p3ty+vWMZHzUPKSimXs97xmkXpw/sWT5GFqMT7c6QyAYVyrJtPy7PZs2ePypcvr/Xr12vv3r0KDQ3VyJEjHYJkly5d9O2332rp0qWaN29elp2P8ePHO4U5SWrZsqUWLFiguLg4+981GzduVGRkpDp06KAPPvhAx48f1x133KH09HRt2bJFQUFBqlOnjv0YoaGhCg0NdTp2gwYN9OSTT+qDDz7QunXr9MADDzht88gjjziFOUmaOXOmJGns2LEO17A1bNhQTz75pCZPnpzl65WVI0eOaPXq1QoMDNTo0aPtYU6S+vbtq3nz5mnr1q25Pl5+rVixQvv27VP37t01cOBAh3UNGjTQiBEjNHbsWC1YsEDPPvusw3o/Pz+NHz/eIcxJ0ldffaWzZ8/queeecwhzktS6dWs98cQTmjZtmlauXKnOnTs7rLedemsLc5LUv39//ec//1FcXJxSUlLsIfCzzz5TcnKynnrqKYcwJ0klSpRQhw4d7D/v3r1bP/30k5o2bapXXnnF4d+q7TEHDBiguXPnavz48bl9+TyaR11DFxMTozlz5mjz5s26cOFCrvZJSUnR+PHjFRMTo+TkZDVt2lQVK1bUqlWr9Pe//12nT58u9DotFovTBbIAAHgz67WrSn9vrJbcCNa0Wt3tyyuV8dEb7UNVtaxv9jtDvXv3zvJ/tj+oN2/eLCnji+WSJUs67d+zZ09J0qZNm5zWValSRY0aNcrycW0hztaJS0xM1P79+9WqVSuHUzKljD/Uk5KS1LJlS6fOZFpamlavXq2JEyfq73//u0aMGKERI0bY9z18+HCWj59VyDt58qT++OMPVa1a1anrJ0ndu3fP8ljZ2bJliySpXbt2WU5U8tBDD+XpePllu57wwQcfzHJ98+bNJWVcp3izBg0aqFq1atkes2PHjnk+ZqtWreydMxsfHx+Fhobqxo0bDn+L//LLL5Kkxx9/PMvHycy27YMPPphlB9s2plld9+mtPKpDV7t2bdWoUUMRERGKiIhwuOA2OwsXLtT+/ftVu3ZtjR492v7NxZIlSzRz5kxNnTpV48aNs29/5coVXblyJcdjli5dOstvsSRp1apVOnjwoAYMGJCHZwYAgOeyXruq9Pdf0+IbVfR5rT+v4apszghzVQJ88z2Rhbe41W0LbF9Q33777Vmuty0/c+aM07rbbrst2+PaQpst0Nmun7v77rvVrFkz+fr6asOGDYqOjs7ydEtJ+uOPP/Tkk09mOQGHjcWS9Sm1WXUmbc81u65ldsuzY3tNsnsd8nq8/LLdeHzw4MEaPHhwttslJiY6LcuuxhMnTkiSU3cuN8fMKiBKsp8am5KSYl92q5k9M7M9zzfeeENvvPFGnmryVh4V6PL6jUtqaqp9ppyBAwc6tKG7dOmi1atXa+/evTp06JDCw8MlSUuXLlVMTEyOx23Tpo2GDh3qtHzLli367LPPNGjQIPvxAADwZtbkjDC3KCVYM2v++UdlsLmUJrS/Q1UCctGZK23OOH3RaEo7XxNW1G51zV5W628+TS+z6tWr67bbbrNfR7dx40aZTCa1aNFCpUuXVqNGjexhzzYL482B7uWXX9aePXvUuXNnDRkyRBEREQoICFCJEiX05Zdf6u9//3u2gd7Pz89pmW3b7J5rXq9bvNXxikt6erqkjE5hTpcEZXUz+axeJymjMypl/N1rmyQlK3fddZfTsry+Hrnd3lZTixYtsjwV1yYoKChPj+/JPCrQ5dW+fftksVhUpUoVhYWFOa1v0aKFjh49qq1bt9oDWK9eveynJWQnq3+w69at00cffaRnnnlGbdu2LZT6AQAwMmvyNaW//7q+SamqLyP+vD6nirmUJrQPVXBAqRz2/pPJZCq0a9E8VdWqVSX92f242cmTJyXJ6SbYudGyZUt98803iouL04YNGxQZGWn/Y7tVq1b64IMPdPToUW3evFkVKlRQZGSkfd+rV69qzZo1qly5sj7++GOn00GPHTuW53psz9XWfbpZdsuzU6VKFUl/vkY3y255YbN1xPr165flqab5PebBgwc1fPjwIp1gpFq1ajp8+LCOHj3qMP7ZbStJUVFRTtcKImteHehs09ZmFeYk2UNc5ultM18Im1srV67U559/riFDhmQ57WtOMk8DbOPr66u3335bkpxmMSpKtoteb579CsbH2HomxtVzecLYpidf04X3XtXXKVU0J6KTfflt5fw1pfedqlLWsaNgtVp14cIF+fj4uLxTUlRszyvzJBO5dat9bF2x7777Tq+88opTcFq4cKF9O9uxbNuYTKYcj3/vvffqm2++0ffff699+/bp6aeftm9/77336oMPPtD06dN16dIlde7c2eG6q6tXryo9PV1VqlRx6iJlPpOqRIkSDjXYXquSJUs61VajRg3ddttt+uOPP7Rt2zan6+gWL16cq+dlY3vtfvrpJ12/ft1pts3vvvsuyxptfzOWLFnSYWxtr+vN22eWVW1t2rTR3LlztXz5cqcJSrJzq8dq06aN1q5dqxUrVujOO+8slGNmNTZt2rTR4cOHNWfOHL355ps5Hr9t27b697//reXLlztN7uIqVqtVJUqUUOXKlR3ef9zlvdijJkXJq4SEBEnKtm1t+3bJtl1+LFmyRNOmTVP//v1Vv359Xbx4URcvXlRSUlK+jwkAgJFZryfr4oSXNOdKkOaE/xnmbi/vr4+yCHMouHvuuUd169bVsWPH9K9//cvhFMalS5cqNjZWZrM52/u95cQ2McqcOXPs18/Z2K6js804fvMM35UqVVK5cuW0b98++8QtUsZpd6+//roOHjyY53ok2e8N99prr+ny5T9nR929e7c+//zzPB0rLCxM9913ny5evKgJEybYT32UpLlz59onTSlqXbt2Va1atTRv3jx9+OGHunHjhsP6lJQUxcbG5ngt4s2eeOIJVaxYUe+//77mzp3rdGqrxWLR119/bb+/X34NGjRI/v7+mjlzppYsWeKwLj093eEG7U2bNtV9992ndevWacyYMU7XT6anp2vVqlVZTuDjrby6Q5ecnCwp+/OKbeeM27bLj++//17p6en67LPP9Nlnn9mXV65cWVOmTLnl/pMmTcpxfUJCQrFdKG779uHcuXPF8ngoPoytZ2JcPZeRx9Z6/brSJ4/XvGvBmhf252ljt5UtpdfbhciUnKRzWXzsWq1WpaenKzU11WM7dLZv+1NTU/O8b272+fDDD9W7d2+9//77Wrp0qerXr6+TJ09qy5Yt8vHx0cSJExUUFGQ/lu1aJqvVmuPx77jjDlWrVk2nTp2SyWRSs2bN7Nv7+vqqUaNG9tDTvHlzp2MNHjxY//rXv9S9e3fdc889CgwM1LZt23Tu3Dn1799fX3zxhX3sbWx/+6SlpWVZ27PPPqsffvhBW7duVfPmzXX33XfLYrFo3bp1euSRRzRz5sxbPq/M3nzzTXXv3l0zZszQ6tWr7TcW37Ztm/r166dZs2Y51WgLfmlpafZ6U1NT7a/rzdtnll1t06ZNU9++fTVhwgR99tlnqlu3rgICAvTHH3/o4MGDunTpkqZPn67atWvbHzunxwoICND06dPVv39/jRgxQu+++67q1KkjPz8/nTx5UvHx8bp69aqWL19uPx33VsfMamyqV6+uiRMnasSIERo4cKAiIyNVp04dXbp0Sfv27dPp06cdTl398MMP9eijj+rTTz/V119/rfr166tixYo6ffq0Dh48qPPnz+u1115TkyZNsh2zwmR7/zl37pzD+0/m92KTyZTtRDFFzasD3a2CUGEEpdyENgAAvIE15brSpkzQvOQq+jrsz/tOhZQtpQkdqiuotFf/WVLk6tatq+XLl+v999/Xzz//rNjYWJUtW1YdO3bUc889l+XEF7nVqlUrffPNN6pbt64qVKjgtG7Lli0KDAzM8jqtYcOGqVq1apo2bZq2bNkif39/NW/eXC+99JJ27dqVr3r8/Pw0b948TZo0SYsWLdIPP/yg22+/XSNHjtSzzz5rv09dboWHh2vx4sV666239Msvv2j58uWKjIzUjBkzFBAQoFmzZuWrzryqWbOmfvjhB82YMUPff/+9Nm/eLKvVqipVqqhFixbq2LGj7rvvvjwds1mzZvrxxx/16aef6scff9S6detUsmRJValSRe3bt1enTp3sAbEgunfvrlq1amnq1Klav369li5dqsDAQNWsWdNpMsHKlStr8eLFmjVrlr777jvt2LFDN27cUHBwsBo0aKAHHnjA6X523sxk9eB5gPv06ZPjjcX/+9//KjY2VlFRUXryySed1h85ckQjR45UWFiY/vWvfxV1ufly6tQpOnQoMMbWMzGunsuIY2u9kaK0yW9o7rVgxVS/37789nKlNKF9dVW4RZizWq06f/68KlasSIcOhsPYGlt27z906NyAbUKR8+fPZ7nedn+L4px4JCcWi8V+HnF+ZqECAMAVrDdSlPbRm5p9rYq+qd7OvvyOcr6a0D5UgXTmACDfvPod1HZzw8OHD2e5/tChQw7buVpsbKxiYmLk5+dXbK19AAAKwnrjhtKmvq1ZV6toUfW/2peHlvfV+PahCvT36j9FAKDAvPpdNDIyUmXKlNGZM2d0+PBhp9sX2GbPady4sSvKcxIVFcU97AAAhpER5t7Sf69W0XehbezLq/8vzJUnzAFAgXn1bQt8fHzUsWNHSdKMGTMcZrNcsmSJ/eaHNWvWdFWJDsxms4KDgzndEgDg9qypN5T2yb/0+dVq+u6OP8NcjcCM0ywJcwBQODzq3TQuLk4LFixwWJaamqpRo0bZf+7Zs6dDx61Hjx7atWuX9u/fr+HDhysyMlIJCQmKj49X2bJlNWTIkGKrHwAAT5AR5t7RDEtVxd5xr315eKCfxrUPVTm/kjnsDQDIC48KdElJSYqPj3dYZrVaHZbdfENvX19fjR07VgsXLtTatWu1ZcsWmc1mtWnTRtHR0W4zIQoAAEZgTU1V2qf/1rQrVfX97ffYl0dU8NO4+0NVljAHAIXKowJd27Zt83WNma+vr6KjoxUdHV34RRUiZrkEALgza2qqUqe9q8+uVNXy2++2L6/5vzAXQJgDgELnUYHO0zHLJQDAXVnT0pQ2fZI+uVxVK0Ja2pfXCvLTa/eHKsCXMAcARYFAZyDMcgkAcEe2MDc1qYpW3tbCvrxORT+NbRcqM2EOAIoMgc5AzGazzGazq8sAAMDOmp6mtBnv6aOkKvrxtub25ZEV/TT2/lCVKUWYA4CiRKADAAD5Yk1PU+rn72tKUhX9XK2pfXndSn56tR1hDgCKA4EOAADkmTU9TalffKgPL1bR6qpN7MvrVfLTGMIcABQbAh0AAMgTa3q6UmdO0fsXq+iXqnfZl9ev5K8x7UJVulQJF1YHAN6Fd1wDsVgsOnv2rM6ePevqUgAAXsqanq7UWR/pvcRg/VLlzzDXsLK/Xr2fMFfcQkJCnP5Xo0YNNW3aVM8995z27t3rkrp69eqlkJAQHT9+3CWPXxzWr1+vkJAQjRgxolCON3HiRIWEhGjevHkOy3N6LRctWqSOHTsqIiJCISEhatHiz0mJjhw5ooEDB6pBgwa6/fbbFRISovXr1xdKrXAvdOgMhNsWAABcyZqerhtfTtV/EitrfZVG9uV3Bvtr9F9D5edDmHOV3r172//78uXL2rlzpxYuXKjY2Fh9+eWXuueee3LYG0a0fft2Pf/88/Lz81ObNm1Urlw5BQUFSZLS09M1aNAg/fbbb2rcuLHCwsJUokQJ7mPsoQh0BsJtCwAArmK1WnVjzieamFhZG4PvtC9vFOyvUYQ5l3vvvfccfr5x44b+9re/acGCBRo7dqxWrlzpmsJQYO+//76uXbumqlWrOixfsWKF0tPTNWHCBD3yyCMO644fP67ffvtNLVq00DfffFOc5cIFePc1ELPZrODgYL5dAQAUK6vVqpQ5n+rdhEraWPnPMHdXFcKcuypVqpT+9re/SZL27t2rS5cuubgi5FdISIhq1qypUqVKOSw/deqUJCk0NNRpn5zWwfPwDgwAALKVEeY+07sJFbWpckP78sZV/PUKYc6tVa5c2f7faWlpDut2796tCRMmqGPHjmrYsKHCwsLUqlUr/fOf/9Tp06ezPebJkyc1atQo3XPPPQoPD1f9+vUVFRWlDz74QNeuXbtlTUlJSerRo4dCQkL06quvymq1avjw4QoJCdGGDRsctl2yZIn9usCbrx/7+OOPFRISoi+++MK+7PDhw5o4caK6du2qv/zlL6pRo4aaNGmiYcOG6eDBg1nWY7vuLCUlRf/5z3/UunVrhYWFacCAAfZtjh8/riFDhqh+/fqqVauWunXrVqCO54YNG9SrVy/VqlVL9evX18CBA3XgwIFst7/5Grp58+Y5XGvXu3dv++tkW9ezZ09J0vz58+3revXq5XDcvXv36rnnnlOTJk0UFhamxo0b64UXXsjyWr3M1/dt27ZNTzzxhOrXr6+QkBDt3r3bvt3mzZs1cOBA3XnnnQoLC1OLFi00ZswYnT9/3umYI0aMsF/Xt3HjRvXu3Vu1a9dWnTp11K9fP/3+++/ZviY//fSTnnjiCfvjNGvWTAMGDMhyXI4fP66RI0eqRYsWCgsLU8OGDfXMM89oz5492R7faHgXBgAAWbJarUr5arreSaikzZUa2Jc3rZoR5nxL8meEO9u5c6ckKSgoyH5tlc2UKVP06aefKi0tTc2aNVO7du1ktVo1c+ZMde7cOctQt3HjRrVv315ffPGFrFarHnjgATVp0kSJiYn617/+pYSEhBzrOXfunHr16qVNmzbppZde0uuvvy6TyaS7775bkpwm7Mgc8G5eZ/u5ZcuW9mVz587VpEmTdOXKFd15553q0KGDypYtqwULFigqKirbP+DT09M1cOBAffTRR6pevboeeOAB+9lQR44cUZcuXfTtt98qKChIHTp0UHp6uvr376/Fixfn+Hyzsnz5ckVHR2vDhg2qV6+e2rRpo71796pLly46cuRIro4RFham3r17q0aNGpKktm3bqnfv3urdu7d9ne0SnRo1atjXZb5sJzY2Vp07d9bChQsVHBysDh06qHLlyvr666/VqVMn7d+/P8vH3rRpkx5++GGdOHFCbdq0UcuWLVWiRMb7wPTp09WjRw+tWLFCNWrUUIcOHeTv768ZM2aoS5cuOnPmTJbHXLFihfr06aOLFy+qTZs2Cg4O1k8//aQePXpkORHguHHj1K9fP61atUoRERHq1KmTQkNDtX79en388ccO227evFkdOnTQ7NmzZTab1aFDB4WFhen7779X165dtW7duly95u6Oa+gAAIATq9Wq6/M+1zvnKurXSnXty5tV9dff24aqlBuFOavVKsuNdFeXkWfmUiVkMpkK/bhJSUnavn27Ro0aJUl6/vnnnbZ57LHH9Nprr6lKlSr2Zenp6Xr//ff17rvv6p133tGkSZPs6y5evKhBgwYpKSlJr732mp5++mmH2jdu3Kjy5ctnW9Px48f1yCOP6OjRo5owYYKeeuop+7pWrVpJklOHbsOGDQoLC9PJkye1YcMGRUdH2+vcsmWLgoKCVKdOHfv2Dz74oPr27WsPOjbz5s3Tiy++qLFjx2r+/PlOtf3xxx/y9fXVmjVrVK1aNYd1r7zyihISEvTkk09qwoQJ9vAyZ84cvfzyy9k+36xcuXJFL730ktLS0jRlyhR1795dkpSamqqXXnopy9qy0rx5czVv3lwjRozQkSNHNHToUHsotq1fv369Vq1apWbNmjldX3ns2DENHz5c/v7+mjt3rkMonj9/vkaMGKEXX3xRsbGxTo89b948jRo1SkOGDHFY/uuvv+q1115TSEiIPv/8c9WrV09Sxu/me++9p3fffVdjxozRp59+6nTMadOm6cMPP7S/Hmlpafq///s/LV26VP/9738dXucFCxbo008/VbVq1TRz5kz740jS1atXFRcXZ//58uXLevbZZ5WcnKxPPvlEXbp0sa9bs2aNnnzySQ0bNkwbNmyQr69vTi+52yPQAQAAB1arVdfn/1f/OhukuIqR9uXNq/lrZJvqKlWy8ENIQVhupOux+fGuLiPPZveupQDfwrkBe0hIiNOySpUqOQSHzO69916nZSVKlNALL7ygL7/8UsuXL3dYN2fOHJ0/f17t27fXM88847Rv5lBws/3796tv375KSEjQBx98oB49ejisDw0NVUhIiOLi4pScnCx/f38lJibq999/11NPPaVdu3Y5hL3du3crKSlJnTt3dgiVTZo0UVaio6M1d+5cbdiwQUlJSSpXrpzTNv/85z+dwtyRI0e0evVqBQYGavTo0fYwJ0l9+/bVvHnztHXr1myf982+++47JSYmqnXr1g5j4uPjo9dee01Lly6VxWLJ9fHya9q0abp27Zpee+01p3Hr3bu3li1bpmXLlmnXrl1q2LChw/rIyEgNHjzY6ZhTpkxRenq63nnnHYeQZTKZNGLECC1btkzff/+9EhMTnbrF3bt3d3g9SpYsqWHDhmnp0qXatGmTw7YffvihpIwuXebHkaQyZco4/Lv+6quvdPbsWT333HMOYU6SWrdurSeeeELTpk3TypUr1blz5+xeLkMg0BmIxWKx/6IzMQoAoChYrVYlx8zU22cqaHvFP7sfLav56yU3DHPIkPm2BSkpKTpx4oS2bdumCRMmqEqVKvYuWGaJiYlasWKF9u3bp6SkJPt1dqmpqbp48aIuXLigChUqSJJ++eUXSdLjjz+ep7ri4uL0yiuvKDk5WdOnT1f79u2z3K5ly5ZasGCB4uLidPfdd2vDhg2yWq1q1aqVAgIC9MEHH+j48eO644477KdbZvWcLBaLVqxYod9++00XL17UjRs3JElnz56V1WrV0aNHnUKKyWRShw4dnI61ZcsWSVK7du1UpkwZp/UPPfRQngKdLZx07drVaV1gYKDatGmjpUuX5vp4+WUbywcffDDL9c2aNdOyZcu0fft2p9eqffv2Tl3l9PR0rV27VgEBAVl+UWAymdSsWTPt3r1bO3fudJqxvU2bNk77hIeHS5LDaZqnT59WfHy8KlSooKioqFs+zzVr1kiSOnbsmOX65s2ba9q0adqxYweBDsWH+9ABAIqS1WpV8jez9dbpCtoRVNu+/O7b/PW3NtXlU4Iw565uPq1Oyuhk9ezZU4899phWrVrlMOPhokWLNHLkyBw7QhaLxR7o/vjjD0lyOp3xVoYNG6bU1FRNnTo12zAnSXfffbcWLFigDRs22AOdyWRSy5Yt7YFu/fr19uvPJOdAt3btWg0ZMiTLCThsrly54rSsUqVK8vPzc1puCxO33XZblsfKqiuaE9vxstsvu8cpbLZJT/7yl7/kuF1iYqLTsqxqv3Dhgv3f0a1m1czqmDd3RqWMmd2ljC8nbPL6b/DEiROS5NSdy01NRkOgMxDuQwcAKCpWq1XJC+fqzVPltTOoln35Pbf560U3D3PmUiU0u3etW2/oZsylivY6xAYNGujxxx/Xxx9/rM8//1xjx46VlPGH7gsvvCCr1apx48bp/vvvV9WqVVW6dGlJUrdu3fTrr7/KarU6HTOv1/w99NBDWrBggd599121bNky2zOMbKf+2cLaxo0bFRkZqaCgIDVr1ky+vr7asGGDevfurS1btqhChQqKjPzzdGCLxaL/+7//04ULFzRixAh1795dt99+u/z9/WUymTR06FAtWrQoy+eUVZiTZN+2sK5zLOzj5Vd6erpMJpPTrJc3y3x9ok1Wr5WtsxsQEKBOnTrleMzbb7/daVlRvR62urp06WL/t52Vu+66q0gevzgR6AzEbDbbv7EAAKCwWK1WXfv2K73xRzntrlDTvvy+kNJ6oXWoSrpxmJMy/iAsrGvRPM0dd9whSQ7T9v/4449KSUnRs88+q6efftppn2PHjjktu+2223TgwAEdPnxYNWvWdFqfnZdffllVq1bVlClTFB0drfnz56tSpUpO29WoUUO33Xab4uLidOrUKe3bt89+64DSpUurUaNG2rBhg3777TddunRJnTp1cggCmzZt0oULF9S5c+csJys5evRormu2sU0Yc/LkySzXZ7c8O7Ybg9s6RzezdaCKWrVq1XTkyBGNHz9eZcuWLfDxgoKC5OfnJx8fnyw7xYXF1sHM7Wyg1apV08GDBzV8+HCn6+08jftMUQUAAFzi2rfz9MYJxzDXOsTfEGEOObOFs8zXgNluMp7VKX4bN27UuXPnnJbfd999kqTZs2fnuYZXXnlFgwcP1u+//67o6OhsT3Fr2bKlrl+/ro8++khWq9Vh5sZWrVrpxIkT+vrrr+0/Z5bTczp8+LDDvdJyq1mzZpIy7nl29epVp/Xffvttno7XvHlzSRn317vZpUuXtHr16jzXmB+269yWLVtWKMfz8fFRq1atdPHiRW3cuLFQjpmVqlWrqlatWrpw4UKurjW0/Zu9eYIfT0SgAwDAi1399itNOFlWuytE2Je1CfHXiNbVCXMGt3v3bnsAa9eunX25bcKJb775xiGonDp1Sv/4xz+yPNajjz6qoKAgrVixQp9//rnTqYubNm1SUlJStrWMHj1agwYN0r59+9SnT58sQ50tpM2ZM0cmk0ktWrTIcl3mn29+Tt9//73DNXSXLl3SSy+9ZJ8cJS/CwsJ077336uLFi3rzzTeVnv7nrTHyOsOllHEqa2BgoFavXq3vvvvOvjwtLU2vv/56scxwKUnPPvus/P399dprr+mHH35wWn/hwgV98cUXubpRvM3zzz+vEiVKaMSIEdq8ebPT+tOnTzvcBD6/hg4dKkkaO3as073yrl69qrVr19p/fvzxx1WxYkV9+OGHmjdvntO/2atXr2r+/PnF1hktSpxyCQCAl7J8O08TTpTTnsBw+7K/3u6v5+8jzBnNiBEj7P9948YNnThxQnFxcUpPT1eHDh0crpd64IEHVKdOHe3YsUP33HOPmjZtquvXr2v9+vWqX7++mjZt6hRWKlSooI8//lgDBgzQ6NGjNW3aNDVs2FDXrl3T77//rmPHjmnjxo1Z3hLAZuzYsUpPT9e0adP0yCOP6Ouvv1ZgYKB9vS2kJScnq169evYJWSTZr6NLTk5WYGCg6tat63DsRo0aqXXr1lqzZo3uu+8+h3vbVahQQQ8++GC+OjVvvfWWunfvrs8//1xr1qxRw4YNdezYMW3btk39+vXL0yR1ZcuW1TvvvKP/+7//0+DBgzVjxgyFhIRo+/btOn/+vHr06KFvvvkmzzXmVXh4uD788EM9//zzeuqppxQREaFatWrJarXqxIkTio+PV0pKih5++OEcrz3LrGXLlnr99dc1duxYPfzww6pbt67CwsJ0/fp1nTx5UvHx8TKbzerfv3+Bau/du7d27Nihzz//XB06dFDTpk1VrVo1nTlzRrt371aDBg3sHcjAwEBNnz5d/fv314svvqhJkyapTp068vPzs9d09epVLV++vNgmpCkqdOgAAPBClsVfa8KJsoQ5DzF//nz7/7777jsdOHBALVq00MSJEzVjxgyHe6j5+vrqm2++0RNPPCE/Pz/9+OOPOnDggAYMGKCvvvpKpUqVyvIx7rnnHv3www967LHHlJqaquXLlysuLk4VK1bUP//5T1WuXPmWdY4bN04DBgzQb7/9pkcffdR+qqSU0RGzzXh4cwfOdh2dlBEesppIY8aMGRo2bJiCgoL0888/a+fOnerWrZsWL16cY9DMSXh4uBYvXqyuXbvq/PnzWr58uaxWq2bMmKFu3brl+XhRUVGaO3euWrRood27d+vnn39WrVq19N133+V5BtGC6Ny5s1asWKHHH39cqamp+vnnn7VhwwZ7kPvvf/+b59fsqaee0pIlS9SjRw9dunRJK1as0K+//iqTyaR+/fppxowZhVL7hAkTNH36dN17773av3+/li5dqmPHjunee+91uuF5s2bN9OOPP9q7kuvWrdPq1at1+fJltW/fXlOnTlXt2rWzeSTjMFmzmu4HhnHq1KksZ2wqCrY36qzOrYexMbaeiXH1XAUdW8vi+Rp/LEB7A8Psy9rd4a/n7nXfMGe1WnX+/HlVrFjR5bMEFhUfn4wTp1JTU11cCQobY2ts2b3/ZH4vNplMWd6CoThwyqWBcGNxAEBBWZbE6PXjAdqXKczdf4e/nruvukp4aFACAE9GoDMQbiwOACgIS+wCvX7MrH3lM4W52wlzAGBkBDoD4cbiAID8urL0G71+1Kz95WvYl3W4w19DCHMAYGgEOgPhxuIAgPy4snShXj9ShjAHAB6IQAcAgAe7/P0ivX6kjH4vX92+7IFQfw2+lzAHAJ6AQAcAgIe6/P23GnektOIzhbkH7/DX/xHmAMBjEOgAAPBAl5f9L8yVC7Uv6xjqr2cJcwDgUQh0AAB4mMvLF+u1w6V1IFOY6/S/MOep93ADAG9FoAMAwINcXr5Yrx3ycwhznUP9NIgwBwAeiUAHAICHSPphiV475K+D5e6wL+tc3V+D7iHMAYCnItABAOABklbE/i/M3W5f1qW6v54mzAGARyvh6gIAAEDBXFqxVGMP+ulg2cxhzo8wBwBegA6dgVgsFlksFklScHCwi6sBALiDSyszwtzhsiH2ZV2r+2ngPTUIcwDgBQh0BhIbG6uYmBj5+flp1qxZri4HAOBiF3/8Xq8dcAxz3ar7aQBhDgC8BqdcGkhUVJQmT56siRMnuroUAICLnVqyUK/F+zqEuYdqEOa81dWrV/Xpp5+qV69eatSokWrUqKF69eqpa9eu+ve//62TJ0+6usRic/z4cYWEhKhXr16Fcrz169crJCREI0aMcFg+ceJEhYSEaN68eU777Ny5U48++qjq1q2rkJAQValSRceOHZOUMVZjxoxR06ZNFRoaqpCQEP62Q4HQoTMQs9kss9ns6jIAAC52aslC/W17io5kCnPdq/up/92EOW/066+/6plnntGZM2dUunRpNW7cWJUrV1ZSUpJ27Nih9957T1OnTtUXX3yh1q1bF2ttvXr10oYNG7Rx40bdcccdt97BA1y5ckVPPfWUzpw5o1atWum2225TyZIl7X/DvfXWW5oxY4Zq1Kihrl27qlSpUqpfv76Lq4aREegAADCQiz/9oFfjfXQ04Db7su41CHPeas+ePerTp4+Sk5M1dOhQjRgxQmXKlLGvT09P17Jly/TGG2/o1KlTLqzU8zz11FN66KGHnOY12L59u06fPq2ePXvqgw8+kCT5+GT8yZ2amqply5bJ399fK1ascBgrIL8IdAAAGMSFn3/Q2PhSOhpQzb6sRw0/PUGY80pWq1XDhg1TcnKy/va3v+nFF1902qZEiRLq3Lmz7r33Xv3xxx8uqNJzBQUFKSgoyGm5LThXr149y/1OnTqlkJAQwhwKDdfQAQBgAIk/r9CrvzuGuZ6EOa+2atUq7d27V9WqVdOwYcNy3LZcuXKKjIx0WGa1WvX111+rR48eqlu3riIiItS+fXt9/PHHunHjhtMxWrRooZCQjNN858yZo/bt2ysiIkJ/+ctfNHLkSF26dMm+re06tg0bNkiSWrZsqZCQEPv/bEaMGKGQkBCtX79eq1atUq9evezXndmOt2nTJo0aNUrt27dXvXr1FBERodatW+vNN990eMyCOn78uIYMGaL69eurVq1a6tatm1auXJnt9jdfQ2d7zrZr7SZNmmR/vsOGDdPDDz+skJAQWa1WnThxIsvXQ5LOnz+v119/Xffdd5/Cw8NVr149Pf7449q4caNTDZmv7zt79qxeeuklNWnSRKGhofrss88cntvIkSPVokULhYWFqWHDhnrmmWe0Z88ep2POmzfPfl3fyZMnNXToUDVs2FARERHq1KmTfvjhh2xfk99//10vvPCCmjdvrrCwMDVq1Eg9evTQtGnTnLa1WCz6z3/+o/vvv18RERGqU6eOevbsqWXLlmV7fGSNDh0AAG4ucdVKvRpfSscDqtqXPV63vHrdVZUwp4xgknrD6uoy8synlKlA4/fjjz9Kkrp06WI/pS+30tPTNXjwYC1ZskRly5ZVo0aNZDabtW3bNo0fP17r1q3Tf//7X5Uo4fzd/4QJEzR9+nQ1atRIbdu21datWzV79mwdOHBACxYskMlkktlsVu/evbVq1SqdO3dOnTt3znEegEWLFmnOnDlq1KiR/vrXv+ro0aP212b8+PHas2eP6tSpo3vuuUfXr1/X7t27NWXKFK1cuVKLFy8u8BwDR44c0UMPPaSEhASFh4erYcOGOnbsmPr3769+/frl6hi253zkyBFt2bJF9erVs18b17x5c9WsWVMhISGaP3++ypQpo6ioKKdjHDhwQNHR0Tp9+rRq1Kihdu3a6cKFC1q3bp1Wr16tDz74QA8//LDTfufPn1fnzp2VlpamZs2a6fr16ypdurQkafPmzXriiSd0+fJl1alTRx06dNDp06f1/fff66efftLMmTN1zz33OB3z+PHj6ty5s/z8/NS8eXOdO3dOv/76qwYOHKgvv/xSbdq0cdh+8eLFGj58uK5fv646deqoadOmunjxovbv36+xY8fq6aeftm977tw59enTR7///ruqVq2q1q1b69q1a/bj//Of/9Rzzz2Xq9cdBDoAANxa4uof9ervPjpu/jPM9atbXv/3QEMlJCS4sDL3kXrDqmULk1xdRp51fLicSvnmP9Dt3r1bktSwYcM87/vxxx9ryZIlat26tSZPnqyKFStKypiBcciQIVqxYoVmzpyp/v37O+37zTffaPHixWrQoIEkKTExUV27dtWmTZu0bt063XvvvQoKCtJ7772nXr166dy5c3r11VdznBRl9uzZ+uijj/TQQw85rXvhhRfUpEkTBQYG2pddv35dY8aM0ezZs/Xpp5/qhRdeyPNrkNkrr7yihIQEPfnkk5owYYI9yM6ZM0cvv/xyro5he87z5s3Tli1b1LFjR/3tb3+T5HgN3fz58+3bZpaWlqZnn31Wp0+f1uuvv64BAwbYQ+3u3bv1yCOPaOTIkbrvvvtUqVIlh31/+uknderUSZMnT5a/v799+eXLl/Xss88qOTlZn3zyibp06WJft2bNGj355JMaNmyYNmzYIF9fX4djzp8/XwMGDNDYsWPt9U+bNk1jx47V+++/7xDoDh06pBEjRig9PV1Tp05Vt27d7OvS09PtXz7YvPjii/r99981ZMgQjRw5UqVKlZIkHT16VH379tU777yjdu3aqV69erl67b0dp1wCAOCmEtf8pFf3O4a53mG++r8HGtKZgy5cuCBJ9jCWW6mpqZo6daoCAgI0ZcoUh/3LlCmjf//73/Lz89OXX36Z5f4vv/yyPcxJGUHmiSeekJRxemR+3H///VmGOdu6zGFOkvz8/DRu3Dj5+Pho+fLl+XpMmyNHjmj16tUKDAzU6NGjHbqSffv2VdOmTQt0/NxasWKF9u3bp+7du2vgwIEOv+MNGjTQiBEjdPXqVS1YsMBpXz8/P40fP94hzEnSV199pbNnz+rZZ591CHOS1Lp1az3xxBM6ffp0lqeWVq9eXa+++qpD97d///4KDAxUXFycUlJS7Ms/++wzJScn6/HHH3cIc1LGdZwdOnSw/7x792799NNPatq0qV555RV7mMv8mGlpaZo7d+6tXjL8Dx06AADc0PnVP+vV/SV1wlzFvqxPDV/1bRVGmIOkjFNN82P37t1KTExUu3btspzUo3LlygoLC9O+fft07do1+6l7Nlnd+iA8PFySdObMmXzV9MADD+S4/tSpU1qxYoUOHDigK1euKD09XZJUqlQpHT58OF+PabNlyxZJUrt27bKcqOShhx7S1q1bC/QYubFmzRpJ0oMPPpjl+ubNm0uSduzY4bSuQYMGqlatmtNy2zE7duyY7TGnTZumHTt2qHPnzg7rWrVq5RC2pIxOY2hoqHbu3KkLFy6oSpWM96dffvlFkvT4449n+/xsbNs++OCDWb6XNWvWTFLGbKHIHQIdAABu5vyanzVmf0mdNP85HXp0WCn1vTvchVW5L59SJnV8uJyry8gzn1IFC+ZBQUE6ePCgzp8/n6f9jh8/LinjNL2bJ+S42cWLF50C3W233ea0ne0atsxdm7zIqY5PPvlEb7/9dr6PfSu2EJrV87pVbYXJNi6DBw/W4MGDs90uMTHRaVl2NZ44cUKSnLpzuTlmVgFRynqsbzWzZ2a25/nGG2/ojTfeyFNNyBqBDgAAN5KwZpVe3ecY5h4J89WjhLlsmUymAl2LZlT169fXli1btGvXLvXs2TPX+9m6W2FhYbc8ndDPz89pWVF0iLN6HCnjpumvv/66ypUrp3/961+6++67VblyZfv2jRs3zndX0MbW6XR159s2Lu3atcvxNNqaNWs6Lcvu9UtLS5OUEehuDuaZ3XXXXU7L8vp65HZ7W00tWrRQaGhotttl1T1G1gh0AAC4iYS1q/Xq/hIOYe7RMF89QphDFu6//3598cUXWrJkiUaPHp3rmS5tnZfIyEiniTncjW0K+5EjR6pPnz4O665du6azZ88W+DFspw2ePHkyy/XZLS9stnHp16/fLU9BzcsxDx48qOHDhxfpBCPVqlXT4cOHdfToUafbY2S1rSRFRUVp4MCBRVaTN2FSFAAA3EDC2jUas7eETpbJHOZKEeaQrb/+9a+qU6eOTp06pQ8++CDHbS9fvqz9+/dLkho1aqRy5cpp/fr1unz5cpHWaLsGKzU1NV/72+4zl9XpkEuWLMn3dYSZ2a7Z+umnn3T16lWn9d9++22BHyM37rvvPkkq1Puw2Y5Z0Iljcvs4s2fPzvW23G+u8BDoDMRisejs2bOF8m0UAMB9nFv3i0bvNemPMpXty/qGldIjd0e4sCq4O5PJpA8++ED+/v6aOHGi3nrrLadAYrVa9cMPP6hTp072SSb8/Pz07LPP6tKlS3rmmWfs11lltmfPnkIJMrbu18GDB/O1v22ylblz5zrc7Pz333/Xm2++WeD6pIxTT++9915dvHhRb775pv3URynjJtvFMSGKlNGxqlmzpr7++mtNmTLF6ebuKSkpWrp0qfbu3ZvrYz7++OOqWLGiPvzwQ82bN88pAF+9elXz58/XH3/8UaDan376afn7+2vWrFmKjY11WHfzbQuaNGmie++9V+vXr9fYsWNlsVictl+9erU2b95coJq8CadcGkhsbKxiYmLk5+enWbNmubocAEAhOLfuF43ZI53KFOYeCyulPoQ55EKDBg301Vdf6ZlnntHkyZM1ffp0NWnSRJUrV1ZSUpJ27typc+fOyd/f36HLNWzYMMXHx2vRokVq3bq1GjRooJCQECUmJurYsWM6duyYHnzwwWxvJZBbDzzwgObPn6/nnntOrVu3VrlyGZPXvPvuu7nav0+fPvrkk0+0YsUKtW7dWo0aNdLFixe1ceNGPfjgg9q+fXuWgTSv3nrrLXXv3l2ff/651qxZY7+x+LZt29SvX79i+bvLx8dH06dPV9++ffXmm29q+vTpqlu3rgICAvTHH3/o4MGDunTpkn15bgQGBmr69Onq37+/XnzxRU2aNEl16tSRn5+fTp48qfj4eF29elXLly/PdlKY3IiIiNDEiRM1YsQIDRo0SJGRkapTp44uXbqkffv26fTp0w6nrk6ePFmPPvqopk2bppiYGNWvX18VK1bU6dOn7RP9vPbaa/aZPZEzAp2BREVFqW3btq4uAwBQSM6uW6tX95h0qsyfNwl+PKyUehPmkAfNmjXTunXrNGvWLK1cuVJ79+7Vxo0bZTabFR4ern79+unRRx91+IO9RIkSmjJlijp37qy5c+dqx44d2rlzp4KCghQSEqLevXs73U8sPzp37qzXXntNc+bM0cqVK3X9+nVJuQ90QUFBio2N1ZtvvqkNGzZoxYoVuuOOO/TSSy9p8ODBuvvuuwtco5TRCVy8eLHeeust/fLLL1q+fLkiIyM1Y8YMBQQEFNsX6TVr1tQPP/ygGTNm6Pvvv9fmzZtltVpVpUoVtWjRQh07drSfsphbzZo1048//qhPP/1UP/74o9atW6eSJUuqSpUqat++vTp16qTatWsXuPbu3burVq1amjp1qtavX6+lS5cqMDBQNWvW1NChQx22rVy5shYvXqxZs2bpu+++044dO3Tjxg0FBwerQYMGeuCBBwrl35+3MFkL4+RjuMypU6cK5fzx3KhcOePb43PnzhXL46H4MLaeiXF1b2fXrdWYvdLp0n+GuX5hPup1t/MMdjfzxrG1Wq06f/68Klas6PLZCIuKbVKT/F5vBvfF2Bpbdu8/md+LTSZTtrd6KGp06AAAKGZn169zCnNPhPmoZy7CHAAAmTEpCgAAxejs+nUas8dKmAMAFAo6dAAAFJOz69dr9B7pTKYw92R4SfVoRZgDAOQPgQ4AgGKQEeasOlO6on1Z//CSerhVLRdWBQAwOk65BACgiJ3dQJgDABQNAh0AAEXo7Ib1Gv0bYQ4AUDQ45RIAgCJyduMGpzD3VHhJdSfMAQAKCR06AACKwJmNGzR6dzphrpCYTCaZTCalp6e7uhQAXiYtLc3+HuSOCHQAABSyMxs3OnXmBhDmCszX11fJycmuLgOAl7l+/bp8fX1dXUa2OOUSAIBClBHm0nXWP8i+bEB4ST1EmCuw0qVL69KlS5IkPz8/lSxZ0sUVFS6r1erw//AcjK0xpaWl6fr160pOTlb58uVdXU62CHQAABSS0xs3aoxTmCtBmCskPj4+Kl++vK5du6ZLly553B/HJUpknDjFaaWeh7E1JpPJJF9fX5UvX14+Pu4bm9y3MgAADOT0/06zPJcpzA0ML6FurWq7sCrP4+Pjo7Jly0ryvG5H5cqVJUnnzp1zcSUobIytMbnrNXM3I9ABAFBApzdu+l+Yq2BfRpgrekb5Yyu3bM/H054XGFsULQIdAAAFkBHm0h3C3NPhJdSVMAcAKAYEOgAA8un0ps0a9Vu6EghzAAAX4bYFAADkw+lNmzVqdxphDgDgUgQ6AADy6NRm5zD3DGEOAOACnHJpIBaLRRaLRZIUHBzs4moAwDud2rxZo3c5nmY5KNykKMIcAMAFCHQGEhsbq5iYGPn5+WnWrFmuLgcAvM6fYS7QviwjzNVxXVEAAK9GoDOQqKgotW3b1tVlAIBX+mPzFo0hzAEA3AyBzkDMZrPMZrOrywAAr/PH5i0avTtN5zOHuTDCHADA9Qh0AADkwB7m/ALty54NN6kzYQ4A4AYIdAAAZOOPzVudwtz/hZvUiTAHAHATBDoAALLwx+atGrU7TYmEOQCAGyPQAQBwkz+22MJcefuyweEmdSTMAQDcDDcWBwAgk5NbtmrUrpvDnAhzAAC3RKADAOB/Tm79VaNvDnNhUsdWkS6sCgCA7BHoAADQ/8LczlTnMHc3YQ4A4L4IdAAAr3dia5xTmBtCmAMAGACTogAAvJotzF3IFOaGhksPcJolAMAA6NABALzWn2GunH0ZYQ4AYCQEOgCAV7o5zJms6RoaZiXMAQAMhVMuAQBe58TWbc5hLtykDnfXdXFlAADkDYEOAOBVTmzdplG7UnWRMAcA8AAEOgCA1zj+63aN3pWqi75lJWWEuecipPatCHMAAGPiGjoAgFc4/ut2jd55wzHMhUvtW9VzcWUAAOQfgQ4A4PGyDXN3E+YAAMZGoAMAeDTCHADAkxHoAAAeizAHAPB0BDoAgEc6Fucc5p4nzAEAPAyBDgDgcY7FbdeYHc5h7n7CHADAwxDoAAAehTAHAPAm3IcOAOAxjv26Q2OyOM2SMAcA8FR06AAAHiEjzKU4hLlh4VbCHADAoxHoAACGl12Ya3d3fRdXBgBA0eKUyyK2ceNGffvttzp9+rRSUlIUFBSke+65R7169ZKPDy8/ABTUsV93aPTOG7pEmAMAeCESRRELCAjQQw89pJCQEPn5+enIkSP69NNPdfXqVQ0YMMDV5QGAof0Z5gIkEeYAAN7HowLdoUOHtHPnTh04cEDx8fG6cOGCSpUqpdmzZ+e4X0pKihYtWqR169YpISFBAQEBatSokaKjo1WxYsUC1dSgQQOHn4ODg7Vnzx7t2rWrQMcFAG93LG4nYQ4A4PU8KtDFxMRo69atedonJSVF48eP1/79+1WhQgU1bdpU586d06pVqxQXF6cJEyaoatWqhVbjiRMntH37dt15552FdkwA8DZH43ZqzI4UwhwAwOt5VKCrXbu2atSooYiICEVERGjQoEG33GfhwoXav3+/ateurdGjR8vf31+StGTJEs2cOVNTp07VuHHj7NtfuXJFV65cyfGYpUuXVvny5R2W9evXT2lpaUpNTVX79u3Vv3//vD9BAECWYW54uFV/JcwBALyQRwW67t2752n71NRULVu2TJI0cOBAe5iTpC5dumj16tXau3evDh06pPDwcEnS0qVLFRMTk+Nx27Rpo6FDhzos+/e//62UlBQdPHhQc+bMUWBgoPr06ZOnegHA2x2N26nRO1OURJgDAECShwW6vNq3b58sFouqVKmisLAwp/UtWrTQ0aNHtXXrVnug69Wrl3r27JnjcU0mk9My22mboaGhMplMmjp1qrp16+YQIgEA2bOHuVIZYa7E/06zJMwBALyZVwe6o0ePSlKWYU6SPcTZtpOkEiUK79Z9aWlpt9zmxRdfdFrm6+urt99+W5JUqVKlQqvnVmy3WahcuXKxPSaKB2PrmTxpXPet3+wU5v5Rv7S6PNDcxZW5hieNLf7EuHouxtYzucu4enWgS0hIkKRsZ7IMCgpy2C4/YmJiVKtWLVWpUkVWq1UHDhzQ7Nmz1bRpU5nN5nwfFwC8xb71mzVi3XnCHAAAWfDqQJecnCxJ8vPzy3K97XRI23b5kZKSohkzZuj8+fMqWbKkKleurKioKHXu3DlX+0+aNCnH9QkJCbJarfmuLy9s3z6cO3euWB4PxYex9UyeMK5H4nZqzE2dueHh6WpxV5ihn1dBecLYwhnj6rkYW8+UeVxNJpOqVavmkjq8OtDdKggVRlDq27ev+vbtW+DjAIC3yS7Mtb27wS32BADAexTeBWEGVLp0aUnS9evXs1xvW87EJQBQvI5s20WYAwAgF7w60NkmFDl//nyW6xMTEx22AwAUvSPbdmnMDsIcAAC54dWnXFavXl2SdPjw4SzXHzp0yGE7V7NYLLJYLJKk4OBgF1cDAIXvzzCXMWkUYQ4AgJx5daCLjIxUmTJldObMGR0+fNjp9gWbNm2SJDVu3NgV5TmJjY1VTEyM/Pz8NGvWLFeXAwCF6vC2XXqVMAcAQJ549SmXPj4+6tixoyRpxowZDrNZLlmyREePHlVkZKRq1qzpqhIdREVFafLkyZo4caKrSwGAQkWYAwAgfzyqQxcXF6cFCxY4LEtNTdWoUaPsP/fs2dOh49ajRw/t2rVL+/fv1/DhwxUZGamEhATFx8erbNmyGjJkSLHVfytms5l71wHwOFmFuRFhaWpzd0MXVwYAgPvzqECXlJSk+Ph4h2VWq9VhWVJSksN6X19fjR07VgsXLtTatWu1ZcsWmc1mtWnTRtHR0UyIAgBFKNswdw9hDgCA3PCoQNe2bVu1bds2z/v5+voqOjpa0dHRhV8UACBLhDkAAArOowKdp2OWSwCe4vC23U5h7oWwNLUmzAEAkCcEOgNhlksAnuDwtt0as+O6LhPmAAAoMAKdgURFReXrlFIAcBeEOQAACheBzkCY5RKAkR3atluv3hTmXgxL032EOQAA8s2r70MHACgeh7YT5gAAKAp06AAARerQ9t16dXvmMJemF8PSCXMAABQCOnQAgCJzKfGiXt2eclOYozMHAEBhIdAZiMVi0dmzZ3X27FlXlwIAuVI+KFC9q6ZKyhzm7nRxVQAAeA5OuTQQblsAwIge6thcpuWbVSHAnzAHAEAhI9AZCLctAGBU3R5s7uoSAADwSAQ6A+G2BQAAAAAy4xo6AAAAADAoAh0AAAAAGBSBDgAAAAAMikAHAAAAAAbFpCgGYrFYZLFYJEnBwcEurgYAAACAqxHoDIT70AEAAADIjEBnINyHDgAAAEBmBDoD4T50AAAAADJjUhQAAAAAMCgCHQAAAAAYFIEOAAAAAAyKQAcAAAAABkWgAwAAAACDYpZLA+HG4gAAAAAyI9AZCDcWBwAAAJAZgc5AuLE4AAAAgMwIdAbCjcUBAAAAZMakKAAAAABgUAQ6AAAAADAoAh0AAAAAGBSBDgAAAAAMikAHAAAAAAZFoAMAAAAAgyLQAQAAAIBBcR86A7FYLLJYLJKk4OBgF1cDAAAAwNUIdAYSGxurmJgY+fn5adasWa4uBwAAAICLEegMJCoqSm3btnV1GQAAAADcBIHOQMxms8xms6vLAAAAAOAmmBQFAAAAAAyKQAcAAAAABkWgAwAAAACDItABAAAAgEER6AAAAADAoAh0AAAAAGBQBDoAAAAAMCgCHQAAAAAYFIEOAAAAAAyKQAcAAAAABkWgAwAAAACD8nF1Acg9i8Uii8UiSQoODnZxNQAAAABcjUBnILGxsYqJiZGfn59mzZrl6nIAAAAAuBiBzkCioqLUtm1bV5cBAAAAwE0Q6AzEbDbLbDa7ugwAAAAAboJJUQAAAADAoAh0AAAAAGBQBDoAAAAAMCgCHQAAAAAYFIEOAAAAAAyKQAcAAAAABkWgAwAAAACDItABAAAAgEER6AAAAADAoAh0AAAAAGBQBDoAAAAAMCgCHQAAAAAYFIEOAAAAAAyKQAcAAAAABkWgAwAAAACDItABAAAAgEER6AAAAADAoHxcXQByz2KxyGKxSJKCg4NdXA0AAAAAVyPQGUhsbKxiYmLk5+enWbNmubocAAAAAC5GoDOQqKgotW3b1tVlAAAAAHATBDoDMZvNMpvNri4DAAAAgJtgUhQAAAAAMCgCHQAAAAAYFIEOAAAAAAyKQAcAAAAABkWgAwAAbslqtbq6BABwewQ6AADglo4eTNHmX67oYmKqq0sBALfFbQsAAIDbSUuzKn5PspKvWXXmjytq2KS0atT0c3VZAOB26NABAAC3c/xwipKvZZxyWaKEVOW2Ui6uCADcE4EOAAC4FVt3ziY03Fely/AnCwBkhXdHAADgVm7uztWs6+/iigDAfRHoAACA26A7BwB5wzskAABwG3TnACBvCHQAAMAt0J0DgLzjXRIAALgFunMAkHcEOgAA4HJ05wAgf3inBAAALkd3DgDyh0AHAABciu4cAOQf75YAAMCl6M4BQP4R6AAAgMvQnQOAguEdEwAAuAzdOQAoGAIdAABwCbpzAFBwvGsCAACXoDsHAAVHoAMAAMWO7hwAFA7eOQEAQLGjOwcAhYNAV0x27dql6OhoDR061NWlAADgUnTnAKDw8O5ZDBITEzVlyhQ1atTI1aUAAOBydOcAoPD4uLqAwnLo0CHt3LlTBw4cUHx8vC5cuKBSpUpp9uzZOe6XkpKiRYsWad26dUpISFBAQIAaNWqk6OhoVaxYscB1paen6/3331fnzp2VnJyskydPFviYAAAYFd05AChcHhPoYmJitHXr1jztk5KSovHjx2v//v2qUKGCmjZtqnPnzmnVqlWKi4vThAkTVLVq1QLVNWfOHPn7+6tr166aP39+gY4FAIDR0Z0DgMLlMYGudu3aqlGjhiIiIhQREaFBgwbdcp+FCxdq//79ql27tkaPHi1//4wPlSVLlmjmzJmaOnWqxo0bZ9/+ypUrunLlSo7HLF26tMqXLy9JiouL09q1a/XOO+/IZDIV4NkBAGB8dOcAoPB5TKDr3r17nrZPTU3VsmXLJEkDBw60hzlJ6tKli1avXq29e/fq0KFDCg8PlyQtXbpUMTExOR63TZs2Gjp0qM6fP6+PPvpIL7zwgsqVK5e3JwMAgAeiOwcAhc9jAl1e7du3TxaLRVWqVFFYWJjT+hYtWujo0aPaunWrPdD16tVLPXv2zPG4tk7cwYMHlZSUpPHjx9vXWa1WWa1WPfLIIxo0aJDatWtXiM8IAAD3RXcOAIqG1wa6o0ePSlKWYU6SPcTZtpOkEiVy/8HTsGFDvfvuuw7LfvjhB23ZskWjRo1SUFBQXksGAMCw6M4BQNHw2kCXkJAgSdnOZGkLXLbt8qp06dIKDQ11WFauXDn5+Pg4Lc/Jiy++6LTM19dXb7/9tiSpUqVK+aovP3x8Mv65VK5cudgeE8WDsfVMjKvnMtrYpqWl66f9f35BWrt+eYVWD3ZhRe7JaOOK3GNsPZO7jKvXnuuQnJxx2oefn1+W623X1Nm2AwAA+fP7niRZrqRKkkqUMOnOxhVcXBEAeA6v7dBZrdYCrc+PPn36qE+fPnnaZ9KkSTmuT0hIKJJas2L79uHcuXPF8ngoPoytZ2JcPZeRxjYtzaptm5PsP4eGl9K15Iu6xvelTow0rsgbxtYzZR5Xk8mkatWquaQOr+3QlS5dWpJ0/fr1LNfblmee/RIAAOQN184BQNEq1EC3fft2xcbG6uzZs4V52CJhu/bs/PnzWa5PTEx02A4AAOQNM1sCQNEr1FMu16xZo3Xr1qlKlSoKDnbvi52rV68uSTp8+HCW6w8dOuSwnTuwWCyyWCyS5PavLwAAdOcAoOgVaqA7duyYJKlmzZq33HbatGlKTExUr1697LcIKE6RkZEqU6aMzpw5o8OHDzvdvmDTpk2SpMaNGxd7bdmJjY1VTEyM/Pz8NGvWLFeXAwBAtujOAUDxKNR31osXL6pkyZIKDAy0L9u2bZtOnTrltG29evX066+/av369YVZQq75+PioY8eOkqQZM2Y4zGa5ZMkSHT16VJGRkbkKp8UlKipKkydP1sSJE11dCgAAOaI7BwDFo1A7dNeuXbNPNmIzdepUJSUl6auvvnJYXq9ePUnS/v37C+Wx4+LitGDBAodlqampGjVqlP3nnj17OnTcevTooV27dmn//v0aPny4IiMjlZCQoPj4eJUtW1ZDhgwplNoKi9lsltlsdnUZAADkiO4cABSfQg105cqVU1JSksOy9PT0LKfVL1++vEqUKFFo07cmJSUpPj7eYZnVanVYdnNtvr6+Gjt2rBYuXKi1a9dqy5YtMpvNatOmjaKjo5kQBQCAfKA7BwDFp1ADXXBwsBITE3XkyBHVqFFDaWlp9kk8Ll++rLJly9q3NZlM8vf31+XLlwvlsdu2bau2bdvmeT9fX19FR0crOjq6UOoAAMCb0Z0DgOJVqO+wLVq0kCQtWrRIkrRhwwalp6dLcp5NMiUlRVevXlWJErzJ55bFYtHZs2cNcVsIAIB3ojsHAMWrUDt0bdq00cKFC7VhwwYdOHBAFy5ckCT95S9/0bfffqs777zTvu3mzZslSUFBQYVZgkdjlksAgDujOwcAxa9QA53ZbNbLL7+sf//73/Zr4yIjI/XMM89oxIgRevvtt/XAAw8oKSlJX375pSSpdu3ahVmCR4uKisrXaaUAABQHunMAUPwKNdBJGQFt8uTJ2rVrl9LT09W4cWP5+PioT58+mj17trZt22bf1mQy6cEHHyzsEjwWs1wCANwV3TkAcI1CD3SS5Ofnp6ZNmzos69atm8xms5YsWaLTp08rMDBQjz76qFvd5w0AAOQP3TkAcI0iCXTZuf/++3X//fcX50MCAIAiRncOAFyHd1sAAFAgdOcAwHWKtUOHgrFYLPb7+gUHB7u4GgAA6M4BgKsR6AyE2xYAANwN3TkAcC0CnYFw2wIAgDuhOwcArkegMxBuWwAAcCd05wDA9fgaDQAA5BndOQBwD7zzAgCAPKM7BwDugUAHAADyhO4cALgP3n0BAECe0J0DAPfBpCgGwn3oAACuRncOANwLgc5AuA8dAMDV6M4BgHsh0BkI96EDALgS3TkAcD8EOgPhPnQAAFeiOwcA7oev1QAAwC3RnQMA98Q7MQAAuCW6cwDgngh0AAAgR3TnAMB98W4MAAByRHcOANwXgQ4AAGSL7hwAuDdmuTQQbiwOAChudOcAwL0R6AyEG4sDAIoT3TkAcH8EOgPhxuIAgOJEdw4A3B+BzkC4sTgAoLjQnQMAY+CdGQAAOKE7BwDGQKADAAAO6M4BgHHw7gwAABzQnQMA4yDQAQAAO7pzAGAsvEMDAAA7unMAYCwEOgAAIInuHAAYEe/SAABAEt05ADAiAh0AAKA7BwAGxY3FDcRischisUiSgoODXVwNAMCT0J0DAGMi0BlIbGysYmJi5Ofnp1mzZrm6HACAh6A7BwDGRaAzkKioKLVt29bVZQAAPAzdOQAwLgKdgZjNZpnNZleXAQDwIHTnAMDYeMcGAMCL0Z0DAGMj0AEA4KXozgGA8fGuDQCAl6I7BwDGR6ADAMAL0Z0DAM/AOzcAAF6I7hwAeAYCHQAAXobuHAB4Dt69AQDwMnTnAMBzEOgAAPAidOcAwLPwDg4AgBehOwcAnoVABwCAl6A7BwCeh3dxAAC8BN05APA8Pq4uALlnsVhksVgkScHBwS6uBgBgJHTnAMAzEegMJDY2VjExMfLz89OsWbNcXQ4AwEDozgGAZyLQGUhUVJTatm3r6jIAAAZDdw4APBeBzkDMZrPMZrOrywAAGAzdOQDwXHw9BwCAB6M7BwCejXd0AAA8GN05APBsBDoAADwU3TkA8Hy8qwMA4KHozgGA5yPQAQDggejOAYB34J0dAAAPRHcOALwDgQ4AAA9Ddw4AvAfv7gAAeBi6cwDgPQh0AAB4ELpzAOBdeIcHAMCD0J0DAO9CoAMAwEPQnQMA78O7PAAAHoLuHAB4HwIdAAAegO4cAHgn3ukBAPAAdOcAwDsR6AAAMDi6cwDgvXi3BwDA4OjOAYD38nF1Acg9i8Uii8UiSQoODnZxNQAAd0B3DgC8G4HOQGJjYxUTEyM/Pz/NmjXL1eUAANwA3TkA8G4EOgOJiopS27ZtXV0GAMBN0J0DABDoDMRsNstsNru6DACAm6A7BwDgazwAAAwoLS2d7hwAgEAHAIAR/b4nie4cAIBABwCA0aSlpWvnrxfsP9OdAwDvxbs/AAAG8/ueJFmupEqiOwcA3o5ABwCAgaSlWenOAQDs+AQAAMBAjh9OoTsHALAj0AEAYBDcdw4AcDM+BQAAMAjH+86Z6M4BALixOAAARnBzd652/XIqXcbkwooAAO6ADh0AAAZwc3fuzsYVXFwRAMAdEOgAAHBzWXXnAsqWcmFFAAB3QaADAMDNOXbnRHcOAGBHoAMAwI1lNbMl3TkAgA2BDgAAN3Zzd46ZLQEAmRHoAABwU9x3DgBwK3wqAADgpujOAQBuhUAHAIAbojsHAMgNPhkAAHBDdOcAALlBoAMAwM3QnQMA5BafDgAAuBm6cwCA3CLQAQDgRujOAQDygk8IAADcCN05AEBeEOgAAHATdOcAAHnFpwQAAG6C7hwAIK8IdAAAuAG6cwCA/PBxdQGe7uuvv1ZMTIzT8smTJys4ONgFFQEA3BHdOQBAfhDoikHFihX11ltvOSwrV66ci6oBALgbunMAgPzyqEB36NAh7dy5UwcOHFB8fLwuXLigUqVKafbs2Tnul5KSokWLFmndunVKSEhQQECAGjVqpOjoaFWsWLHAdZUoUUKBgYEFPg4AwDPRnQMA5JdHBbqYmBht3bo1T/ukpKRo/Pjx2r9/vypUqKCmTZvq3LlzWrVqleLi4jRhwgRVrVq1QHVdvHhRgwcPltVqVWhoqHr27Kk6deoU6JgAAM9Adw4AUBAeFehq166tGjVqKCIiQhERERo0aNAt91m4cKH279+v2rVra/To0fL3z/hWdMmSJZo5c6amTp2qcePG2be/cuWKrly5kuMxS5curfLly0uSatWqpaFDhyokJERXr17VypUr9eqrr2rUqFG68847C/BsAQCegO4cAKAgPCrQde/ePU/bp6amatmyZZKkgQMH2sOcJHXp0kWrV6/W3r17dejQIYWHh0uSli5dmuUkJ5m1adNGQ4cOlSTdddddDuvq1q2rhIQEfffddwQ6APBydOcAAAXlUYEur/bt2yeLxaIqVaooLCzMaX2LFi109OhRbd261R7oevXqpZ49e+Z4XJPJlOP6iIgI/frrr/kvHADgEejOAQAKyqsD3dGjRyUpyzAnyR7ibNtJGROcFNSRI0dyPdnKiy++6LTM19dXb7/9tiSpUqVKBa4nt3x8Mv65VK5cudgeE8WDsfVMjKt7S0tL10/7//x8qV2vvEKr5+52NoytZ2JcPRdj65ncZVy9OtAlJCRIUrbhKigoyGG7/Jg5c6YaN26s4OBgXb16VStWrNBvv/2ml19+Od/HBAAY3+97kmS5kipJKlHCpDubVHBxRQAAI/LqQJecnHHdgp+fX5brbdfU2bbLj8TERE2ePFlJSUkqU6aMQkNDNWbMGDVo0CBX+0+aNCnH9QkJCbJarfmuLy9s3z6cO3euWB4PxYex9UyMq/tKS7Nq2+Yk+8+h4aV0LfmiruXy44ax9UyMq+dibD1T5nE1mUyqVq2aS+rw6kB3qyBUGEFpxIgRBT4GAMCzcO0cAKCwePVUWqVLl5YkXb9+Pcv1tuWZZ78EAKAgmNkSAFCYvPoTxDahyPnz57Ncn5iY6LAdAAAFRXcOAFCYvPqUy+rVq0uSDh8+nOX6Q4cOOWznahaLRRaLRZIUHJy7mdAAAO6D7hwAoLB5daCLjIxUmTJldObMGR0+fNjp9gWbNm2SJDVu3NgV5TmJjY1VTEyM/Pz8NGvWLFeXAwDII7pzAIDC5tVfC/r4+Khjx46SpBkzZjjMZrlkyRIdPXpUkZGRqlmzpqtKdBAVFaXJkydr4sSJri4FAJBHdOcAAEXBozp0cXFxWrBggcOy1NRUjRo1yv5zz549HTpuPXr00K5du7R//34NHz5ckZGRSkhIUHx8vMqWLashQ4YUW/23YjabZTabXV0GACAf6M4BAIqCRwW6pKQkxcfHOyyzWq0Oy5KSkhzW+/r6auzYsVq4cKHWrl2rLVu2yGw2q02bNoqOjmZCFABAgdGdAwAUFY8KdG3btlXbtm3zvJ+vr6+io6MVHR1d+EUBALwe3TkAQFHxqEDn6ZjlEgCMh+4cAKAoEegMhFkuAcB46M4BAIoSgc5AoqKi8nVKKQDANejOAQCKGoHOQJjlEgCMhe4cAKCo8TUhAABFgO4cAKA48MkCAEARoDsHACgOnHJpIMxyCQDGQHcOAFBcCHQGwiyXAGAMdOcAAMWFQGcgzHIJAO6P7hwAoDgR6AyEWS4BwP3RnQMAFCe+MgQAoJDQnQMAFDc+ZQAAKCR05wAAxY1ABwBAIaA7BwBwBT5pAAAoBHTnAACuwKQoBsJ96ADAPdGdAwC4CoHOQLgPHQC4J7pzAABXIdAZCPehAwD3Q3cOAOBKBDoD4T50AOB+6M4BAFyJrxABAMgnunMAAFfjUwcAgHyiOwcAcDUCHQAA+UB3DgDgDvjkAQAgH+jOAQDcAYEOAIA8ojsHAHAXfPoAAJBHdOcAAO6C2xYYiMVikcVikSQFBwe7uBoA8E505wAA7oRAZyCxsbGKiYmRn5+fZs2a5epyAMAr0Z0DALgTAp2BREVFqW3btq4uAwC8Ft05AIC7IdAZiNlsltlsdnUZAOC16M4BANwNXysCAJALdOcAAO6ITyIAAHKB7hwAwB0R6AAAuAW6cwAAd8WnEQAAt0B3DgDgrgh0AADkgO4cAMCd8YkEAEAO6M4BANwZgQ4AgGzQnQMAuDvuQ2cgFotFFotFkhQcHOziagDA89GdAwC4OwKdgcTGxiomJkZ+fn6aNWuWq8sBAI9Gdw4AYAQEOgOJiopS27ZtXV0GAHgFunMAACMg0BmI2WyW2Wx2dRkA4PHozgEAjIJPJwAAbkJ3DgBgFAQ6AAAyoTsHADASPqEAAMiE7hwAwEgIdAAA/A/dOQCA0fApBQDA/9CdAwAYDYEOAADRnQMAGBOfVAAAiO4cAMCYCHQAAK9Hdw4AYFR8WgEAvB7dOQCAURHoAABeje4cAMDI+MQCAHg1unMAACPzcXUByD2LxSKLxSJJCg4OdnE1AGB8dOcAAEZHoDOQ2NhYxcTEyM/PT7NmzXJ1OQBgeHTnAABGR6AzkKioKLVt29bVZQCAR6A7BwDwBAQ6AzGbzTKbza4uAwA8At05AIAn4KtIAIDXoTsHAPAUfHoBALwO3TkAgKcg0AEAvArdOQCAJ+ETDADgVejOAQA8CYEOAOA16M4BADwNn2IAAK9x/BDdOQCAZyHQAQC8QlqaVfF76c4BADwLn2QAAK9Adw4A4IkIdAAAj0d3DgDgqfg0AwB4PLpzAABPRaADAHg0unMAAE/GJxoAwKPRnQMAeDICHQDAY9GdAwB4Oj7VAAAei+4cAMDTEegAAB6J7hwAwBvwyQYA8Eh05wAA3oBABwDwOHTnAADewsfVBSD3LBaLLBaLJCk4ONjF1QCA+6I7BwDwFgQ6A4mNjVVMTIz8/Pw0a9YsV5cDAG6J7hwAwJsQ6AwkKipKbdu2dXUZAODW6M4BALwJgc5AzGazzGazq8sAALdFdw4A4G34lAMAeAy6cwAAb0OgAwB4BLpzAABvxCcdAMAj0J0DAHgjAh0AwPDozgEAvBWfdgAAw6M7BwDwVgQ6AICh0Z0DAHgzPvEAAIZGdw4A4M0IdAAAw6I7BwDwdnzqAQAMi+4cAMDbEegAAIZEdw4AAAIdAMCg6M4BAECgAwAYEN05AAAy8OkHADAcunMAAGQg0AEADIXuHAAAf+ITEABgKHTnAAD4E4EOAGAYdOcAAHDEpyAAwDDozgEA4IhABwAwBLpzAAA445MQAGAIdOcAAHBGoAMAuD26cwAAZI1PQwCA26M7BwBA1gh0AAC3RncOAIDs8YkIAHBrdOcAAMgegQ4A4LbozgEAkDMfVxfg6S5fvqyvvvpKW7du1eXLl1WhQgU99NBDeuCBB1xdGgC4PbpzAADkjEBXhJKTk/Xqq68qKChIw4cPV6VKlXTx4kWlpqa6ujQAcHt05wAAuDWPCXSHDh3Szp07deDAAcXHx+vChQsqVaqUZs+eneN+KSkpWrRokdatW6eEhAQFBASoUaNGio6OVsWKFQtU03fffaeUlBT94x//UKlSpSRJwcHBBTomAHgLunMAANyaxwS6mJgYbd26NU/7pKSkaPz48dq/f78qVKigpk2b6ty5c1q1apXi4uI0YcIEVa1aNd81bdq0SXXq1NF///tfbd68WaVLl1bjxo0VHR0tf3/+MAGA7NCdAwAgdzwm0NWuXVs1atRQRESEIiIiNGjQoFvus3DhQu3fv1+1a9fW6NGj7SFryZIlmjlzpqZOnapx48bZt79y5YquXLmS4zFLly6t8uXLS5JOnz6t06dPq1WrVvr73/+uCxcuaPr06UpMTNQLL7xQgGcLAJ6N7hwAALnjMYGue/fuedo+NTVVy5YtkyQNHDjQoWPWpUsXrV69Wnv37tWhQ4cUHh4uSVq6dKliYmJyPG6bNm00dOhQSZLValXZsmU1ePBglSxZ0v64kyZN0oABA+zBDwDwJ7pzAADknscEurzat2+fLBaLqlSporCwMKf1LVq00NGjR7V161Z7oOvVq5d69uyZ43FNJpP9vytUqKDKlSvbw5wk3X777ZKkc+fOEegAIAt05wAAyD2vDXRHjx6VpCzDnCR7iLNtJ0klSuTtG+LIyEjt2bNH6enp9n3/+OMPSUyOAgBZoTsHAEDeeG2gS0hIkKRsZ7IMCgpy2C4/unbtqg0bNmj69OmKiopSYmKiZs2apXvvvVflypXL1TFefPFFp2W+vr56++23JUmVKlXKd3155eOT8c+lcuXKxfaYKB6MrWcy4rju3XUxU3fOpOb33KaAsqVcXJX7MeLY4tYYV8/F2HomdxlXrw10yckZ3wD7+fllud52TZ1tu/yoUaOG/vnPf2rOnDl6+eWXFRgYqBYtWqhPnz75PiYAeKq0tHTt2HrB/nPt+uUIcwAA3ILXBjqr1Vqg9bnVsGFDvfXWW/nef9KkSTmuT0hIKLRab8X27cO5c+eK5fFQfBhbz2S0cT0Sf11XLamSMq6du72G1TC1FzejjS1yh3H1XIytZ8o8riaTSdWqVXNJHV57YULp0qUlSdevX89yvW0594sDgKLHtXMAAOSP135a2q49O3/+fJbrExMTHbYDABQdZrYEACB/vPaUy+rVq0uSDh8+nOX6Q4cOOWznDiwWiywWiyRmyQTgOejOAQCQf14b6CIjI1WmTBmdOXNGhw8fdrp9waZNmyRJjRs3dkV5WYqNjVVMTIz8/Pw0a9YsV5cDAIWC7hwAAPnntV+B+vj4qGPHjpKkGTNmOMxmuWTJEh09elSRkZGqWbOmq0p0EhUVpcmTJ2vixImuLgUACgXdOQAACsZjOnRxcXFasGCBw7LU1FSNGjXK/nPPnj0dOm49evTQrl27tH//fg0fPlyRkZFKSEhQfHy8ypYtqyFDhhRb/blhNptlNptdXQYAFBq6cwAAFIzHBLqkpCTFx8c7LLNarQ7LkpKSHNb7+vpq7NixWrhwodauXastW7bIbDarTZs2io6OZkIUAChCdOcAACg4jwl0bdu2Vdu2bfO8n6+vr6KjoxUdHV34RQEAskV3DgCAgvOYQOcNmOUSgKegOwcAQOEg0BkIs1wC8BR05wAAKBwEOgOJiorK12mlAOBO6M4BAFB4CHQGwiyXADwB3TkAAAoPX4kCAIoN3TkAAAoXn6IAgGJDdw4AgMJFoAMAFAu6cwAAFD6uoTMQblsAwMjozgEAUPgIdAbCbQsAGBXdOQAAigaBzkC4bQEAo6I7BwBA0SDQGQi3LQBgRHTnAAAoOnyiAgCKFN05AACKDoEOAFBk6M4BAFC0+FQFABSZEiWkhk3KqFxgSbpzAAAUAa6hAwAUGZPJpKohpVTlNh9dvpROdw4AgEJGoDMQ7kMHwKhMJpPKBZZ0dRkAAHgcAp2BcB86AAAAAJkR6AyE+9ABAAAAyIxAZyDchw4AAABAZlydDgAAAAAGRaADAAAAAIMi0AEAAACAQRHoAAAAAMCgCHQAAAAAYFDMcmkg3FgcAAAAQGYEOgPhxuIAAAAAMiPQGQg3FgcAAACQGYHOQLixOAAAAIDMmBQFAAAAAAyKQAcAAAAABkWgAwAAAACDItABAAAAgEER6AAAAADAoAh0AAAAAGBQBDoAAAAAMCjuQ2dwJpPJKx4TxYOx9UyMq+dibD0T4+q5GFvPZDKZXDq2JqvVanXZoyNPLBaLLBaLJCk4ONjF1QAAAABwNTp0BhIbG6uYmBiVK1dO06ZNc3U5AAAAAFyMa+gMJCoqSpMnT9abb76ps2fP2rt1Z8+eVb9+/XT27Nks98tpfXbrslr+j3/8Q//4xz8K8RkVzK2ed3EfMy/75nbb/IxdTuuyW+5OY2vkcc3t9vzOGm9s3WlcJfcaW08f19xsx+9s0R/TyO/FknuNrTuNa173Lerf2byOubuMKx06AzGbzTKbzVmuu379eo775rQ+u3U3L09JSblFhcXvVs+7uI+Zl31zu21+xi6ndVktd7exNfK45nZ7fmfd45iF/TtbHOMqud/Yevq45mY7fmeL/phGfS+W3G9s3Wlc87pvUf/O5mXM3WVc6dABAAAAgEER6DyA2WxWr169su3e5bQ+u3W3OqY7KIoaC3LMvOyb223zM3Y5rWNci3Zcc7s9v7PGG1vGNXuePq652Y6xLfpj8l5ceNxpXPO6b1H/zub37y5XY5ZL5NqLL74oSZo0aZKLK0FhY2w9E+PquRhbz8S4ei7G1jO5y7jSoQMAAAAAg6JDBwAAAAAGRYcOAAAAAAyKQAcAAAAABkWgAwAAAACDItABAAAAgEER6AAAAADAoAh0AAAAAGBQBDoAAAAAMCgCHQAAAAAYFIEOAAAAAAzKx9UFwLt8/fXXiomJcVo+efJkBQcHu6AiFLZdu3ZpwoQJqlSpkqZMmeLqclBAGzdu1LfffqvTp08rJSVFQUFBuueee9SrVy/5+PARYlQ///yz1qxZo2PHjunGjRuqVq2aunTpovvuu8/VpaGA9uzZoyVLlujIkSNKSEhQr1691KdPH1eXhTzYsWOH5syZoxMnTiggIEB//etf1adPH5UoQR/GyIryd5NPYxS7ihUr6q233nJYVq5cORdVg8KUmJioKVOmqFGjRjp58qSry0EhCAgI0EMPPaSQkBD5+fnpyJEj+vTTT3X16lUNGDDA1eUhn3bt2qUmTZroscceU0BAgDZv3qzJkyerZMmSuvvuu11dHgogOTlZt99+u+6991598cUXri4HeXT48GG9/fbbevDBB/X888/rxIkT+vjjj5WWlqbHHnvM1eWhAIryd5NA58EOHTqknTt36sCBA4qPj9eFCxdUqlQpzZ49O8f9UlJStGjRIq1bt04JCQkKCAhQo0aNFB0drYoVKxa4rhIlSigwMLDAx/FW7jqu6enpev/999W5c2clJycT6PLBHce2QYMGDj8HBwdrz5492rVrV4GO603ccVyHDRvm8HO3bt3022+/af369QS6PHDHsW3cuLEaN24sSbesA3lTHOO9ePFihYaGqn///pKk22+/XYmJiZozZ4569uwpf3//onp6Xq04xrYofzcJdB4sJiZGW7duzdM+KSkpGj9+vPbv368KFSqoadOmOnfunFatWqW4uDhNmDBBVatWLVBdFy9e1ODBg2W1WhUaGqqePXuqTp06BTqmN3HXcZ0zZ478/f3VtWtXzZ8/v0DH8lbuOraZnThxQtu3b9edd95ZaMf0dEYYV0myWCwKDQ0t1GN6OqOMLQpHcYz3/v37nU59bty4sb744gsdOnRI9erVK5TnAkdG/10m0Hmw2rVrq0aNGoqIiFBERIQGDRp0y30WLlyo/fv3q3bt2ho9erT9m6AlS5Zo5syZmjp1qsaNG2ff/sqVK7py5UqOxyxdurTKly8vSapVq5aGDh2qkJAQXb16VStXrtSrr76qUaNG8QdiLrnjuMbFxWnt2rV65513ZDKZCvDsvJs7jq1Nv379lJaWptTUVLVv397+7TFuzZ3H1WbVqlU6ePAgp9HmkRHGFoWnOMb7woULqlChgsMxbGc1JSYmFt6TgYPiGNuiRKDzYN27d8/T9qmpqVq2bJkkaeDAgQ5t/S5dumj16tXau3evDh06pPDwcEnS0qVLs5zkJLM2bdpo6NChkqS77rrLYV3dunWVkJCg7777jkCXS+42rufPn9dHH32kF154gWshC8jdxjazf//730pJSdHBgwc1Z84cBQYGMtFCLrnzuErSli1b9Nlnn2nQoEH24yF33H1sUbiKY7yzYvuilC9Mi46rxrawEOhgt2/fPlksFlWpUkVhYWFO61u0aKGjR49q69at9n+cvXr1Us+ePXM87q3egCIiIvTrr7/mv3DkqKjH9eDBg0pKStL48ePt66xWq6xWqx555BENGjRI7dq1K8RnBJvi/J21nTYSGhoqk8mkqVOnqlu3blzPUQSKc1zXrVunjz76SM8884zatm1bKPUje676nIVr5Ge8K1SooAsXLjhsZ/v55s4dXCc/Y1uUCHSwO3r0qCRl+Q9Tkv0fpG07SYUyhe6RI0cKZVIOZK2ox7Vhw4Z69913HZb98MMP2rJli0aNGqWgoKC8loxcctXvrE1aWlqhHQt/Kq5xXblypT7//HMNGTJE99xzTz4qRV65+ncWxSs/412nTh1t375djzzyiH3Ztm3b5OvrSwfdjeRnbIsSgQ52CQkJkpRtuLL9YW7bLj9mzpypxo0bKzg4WFevXtWKFSv022+/6eWXX873MZGzoh7X0qVLO02kUK5cOfn4+DDBQhErjt/ZmJgY1apVS1WqVJHVatWBAwc0e/ZsNW3aVGazOd/HRfaKY1yXLFmiL7/8UgMHDlT9+vV18eJFSRnhgVOni05xjG1ycrJOnz4tKeO0sIsXL+rIkSPy8fHR7bffnu/jIu/yM95dunTRqFGjNHPmTLVr104nT57U119/rU6dOnFGhBvJz9gW5e8mgQ52ycnJkiQ/P78s19veSGzb5UdiYqImT56spKQklSlTRqGhoRozZozT1OgoPMUxrnCN4hjblJQUzZgxQ+fPn1fJkiVVuXJlRUVFqXPnzvk+JnJWHOP6/fffKz09XZ999pk+++wz+/LKlStrypQp+T4uclYcY3vw4EGHiRhWrlyplStXMrYukJ/xDg8P19///nfNnTtXy5YtU9myZfXAAw8oOjq66AtGruVnbIvyd5NABzur1Vqg9bkxYsSIAh8DeVMc43qzPn36MGFGMSiOse3bt6/69u1b4OMg94pjXPnD3jWKY2zr16+vr7/+usDHQcHld7z/8pe/6C9/+UsRVITCkp+xLcrfTU7Mhl3p0qUlSdevX89yvW05LX9jYVw9F2PrmRhXz8XYehfG23O529gS6GBXqVIlSdL58+ezXG+7/4ltOxgD4+q5GFvPxLh6LsbWuzDensvdxpZAB7vq1atLkg4fPpzl+kOHDjlsB2NgXD0XY+uZGFfPxdh6F8bbc7nb2BLoYBcZGakyZcrozJkzWf4D3bRpkySpcePGxV0aCoBx9VyMrWdiXD0XY+tdGG/P5W5jS6CDnY+Pjzp27ChJmjFjhsPMPEuWLNHRo0cVGRmpmjVruqpE5APj6rkYW8/EuHouxta7MN6ey93G1mQtiinu4Bbi4uK0YMEC+8/x8fEymUwO/7h69uzp8O1BSkqKxo0bp/j4eFWoUEGRkZFKSEhQfHy8ypYtqzfeeENVq1Yt1ucBR4yr52JsPRPj6rkYW+/CeHsuo48tty3wYElJSYqPj3dYZrVaHZYlJSU5rPf19dXYsWO1cOFCrV27Vlu2bJHZbFabNm0UHR3NhbtugHH1XIytZ2JcPRdj610Yb89l9LGlQwcAAAAABsU1dAAAAABgUAQ6AAAAADAoAh0AAAAAGBSBDgAAAAAMikAHAAAAAAZFoAMAAAAAgyLQAQAAAIBBEegAAAAAwKAIdAAAAABgUAQ6AAAAADAoAh0AAAAAGBSBDgAAAAAMikAHAAAAAAZFoAMAAAAAgyLQAQAAAIBBEegAAAAAwKAIdAAAAABgUAQ6AACK0aJFi9SnTx89+uijOnDgQJbbxMXFKTo6Wn369NHatWuLuUIAgJEQ6AAAKEYPPfSQGjZsqLS0NL3//vu6du2aw/oLFy7oo48+ktVqVevWrXXvvfe6qFIAgBEQ6AAAKEYmk0nPP/+8ypcvrzNnzuizzz6zr7NarZo8ebKSkpJUtWpVPf300y6sFABgBAQ6AACKWWBgoIYMGSKTyaS1a9dq1apVkqRvv/1Wu3btUsmSJTV8+HD5+/u7tlAAgNsj0AEA4AJ33XWXoqKiJEkzZszQmjVrNG/ePEnSo48+qoiICFeWBwAwCJPVarW6uggAALxRamqqRo8erUOHDtmXNWrUSK+88opMJpMLKwMAGAUdOgAAXMTHx0dDhgyx/1ymTBkNHTqUMAcAyDUCHQAALrRy5Ur7f1+7dk1HjhxxXTEAAMMh0AEA4CK//vqrli1bJkmqXr26rFarpkyZoosXL7q2MACAYRDoAABwAdv95iSpbdu2GjdunCpXrqxLly5pypQp4hJ3AEBuEOgAAChm6enpmjx5si5fvqxq1appwIABKlOmjIYPH66SJUtqx44dWrJkiavLBAAYAIEOAIBi9t1332V5v7natWurV69ekqS5c+c6zH4JAEBWCHQAABSjAwcOONxvLjw83GH9ww8/rPr16ys1NVXvv/++kpOTXVEmAMAgCHQAABSTa9eu6f3331daWpruvPNOde3a1WmbEiVK6LnnnlNAQIBOnTqlGTNmuKBSAIBRcGNxAAAAADAoOnQAAAAAYFAEOgAAAAAwKAIdAAAAABgUgQ4AAAAADIpABwAAAAAGRaADAAAAAIMi0AEAAACAQRHoAAAAAMCgCHQAAAAAYFAEOgAAAAAwKAIdAAAAABgUgQ4AAAAADIpABwAAAAAGRaADAAAAAIMi0AEAAACAQRHoAAAAAMCgCHQAAAAAYFAEOgAAAAAwqP8HY/DylRPyor0AAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Exact difference\n",
"df = np.cos(x)\n",
"\n",
"# Reducing dx and confirm errors\n",
"dxs = [0.01*np.pi, 0.005*np.pi, 0.0025*np.pi, 0.00125*np.pi, 0.000625*np.pi]\n",
"\n",
"# Fowrad difference\n",
"err_fw = []\n",
"for dx in dxs:\n",
" fdf = (f(x+dx) - f(x)) / dx\n",
" err_fw.append(abs((fdf -df)/df))\n",
" \n",
"# Backward difference\n",
"err_bw = []\n",
"for dx in dxs:\n",
" bdf = (f(x) - f(x-dx)) / dx\n",
" err_bw.append(abs((bdf -df)/df)) \n",
" \n",
"# Central difference\n",
"err_ct = []\n",
"for dx in dxs:\n",
" cdf = (f(x+dx) - f(x-dx)) / (2*dx)\n",
" err_ct.append(abs((cdf -df)/df)) \n",
" \n",
"# Plot log-log graph \n",
"plt.loglog(dxs, err_fw)\n",
"plt.loglog(dxs, err_bw)\n",
"plt.loglog(dxs, err_ct)\n",
"\n",
"# Same scale of x and y\n",
"plt.axis('equal')\n",
"\n",
"# Legend\n",
"plt.legend([\n",
" 'Forward difference',\n",
" 'Backward difference',\n",
" 'Central difference',\n",
"])\n",
"\n",
"plt.xlabel('x')\n",
"plt.ylabel('$\\epsilon$')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}